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Abstract. The propagation of extragalactic jets is studied by a series of twelve axisymmetric hydrodynamic 
simulations. Motivated by observational constraints, but unlike most previous simulations, the regime of jet to 
external medium density (r^) from 10“® to 10“^ is explored, for Mach numbers (M) between 2.6 and 26. The 
computational domain contained the bow shocks for the whole simulation time. The bow shocks are found to be 
spherical at source sizes below a critical value ri (blastwave phase), which can reach up to 10 jet radii. After that, 
their aspect ratio rises slowly, as long as the bow shock stays supersonic. The cocoons expand typically to almost 
the same size as the bow shock, unless the Mach number is below approximately three. Low values for the aspect 
ratio and the cocoon-to-bow-shock width ratio is demanded by recent Chandra X-ray observations of the bow 
shock in the archetypical radio galaxy Cygnus A. Therefore, rj < 10“® and M < 6, in this source. The numerical 
work is complemented by an analytic approach for the spherical phase. Extending previous work, the radial force 
balance could be integrated for arbitrary background density and energy input, which results in a global solution. 

The analytic results are shown to be consistent with the numerical work, and a lower limit to ri can be calculated, 
which falls below the numerical results by a few jet radii. It is shown explicitely how a King density distribntion 
changes the discussed aspects of the bow shock propagation. Because the jet head propagates very fast in the 
blastwave phase, it turns out that it is not possible to “frustrate” a jet by a high density environment. This is very 
important for the class of small radio galaxies (compact symmetric objects / GHz peaked sources): They have 
to be young. During its blastwave phase, a powerful jet can transfer typically 10®° erg to the environmental gas. 

This is enough to balance the radiative losses in a cooling flow, if one of the cluster galaxies harbours a powerful 
jet every 10® years. 
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1. Introduction 

Extragalactic jets have been known in the optica l since 
the f amous observation of the jet beam in M 87 ( Curtis 
1918|). The understanding of those collimated flows im- 


prove d considerably with the ad vent of radio astronomy 
(e.g. Baade &: Minkowski 1954 ). During that era, plenty 
of radio jets have been found in extragalactic objects. 
Surrounding the beams, a large cocoon with lower sur¬ 
face brightness was found, which can extend up to tens 
of jet radii (e.g. in Cygnus A [Carilli fc Barthel 19961 ). 
They have been imaged down to the very center of the 
mainly elliptical host galaxies, where a supermassive black 
hole is believed to power the bipolar, relativistic, outflows 
by an accretion process. Jets can propagate up to several 
Mpc (Schoenmakers et al, 2001). According to standard 
jet theory, they displace the medium they propagate into 
and form a bow shock around the system. Although, the 
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presence of bow shocks could be anticipated from some 
ROSAT images ( |Carilli et al. 1994), the Chandra X-ray 
observatory provided the firs t clear view on tho se features 
with satisfactory resolution (^mith et al. 2002). The bow 
shock in Cygnus A is a very fine example. The source is 
believ ed to be located nearly perpendicular to the line of 
sight (Krichbaum et al. 1998). Therefore, one can easily 
measure the true aspect ratio (bow shock diameter in jet 
direction over perpendicular bow shock diameter) from 
the Chandra image (Smith et al. 2002), which turns out 
to be 1.2. Another parameter that can be easily measured 
with the new Chandra data is the relative extention of ra¬ 
dio cocoon and bow shock. Defining this ratio to be taken 
perpendicular to the jet axis through the core renders it 
independent from the viewing angle. One can in principle 
relate these observables to jet parameters via numerical 
simulation results. However, this is only possible if the 
bow shock stays within the computational domain during 
the simulation. Up to now only a few hydrodynamical sim¬ 
ulations in the literature meet this requirement. Two ex- 
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ample s can be found in Reynolds et al. ( 2002 ) and Krause 
(2002|). In the latter case, a jet with a density contrast (jet 


over surrounding medium) 77 = 10“^ was simulated. This 
resulted in an aspect ratio of slightly more than unity, 
similar to Cygnus A. 

Jet simulations have been performed so far mainly 


Camenzind 


with 10 ^ < 77 < 

10 (e.g. Norman et al. 

19821, 1983; 

Clarke et al. 198f 

3; [Kossl & Muller 1 

988; 

jind et al. 

1981; 

Cioffi & Blondin 1992; Clarke 1992 

; Massaglia et al. 

1996; 

Komissarov 

1999; Fregillis et al. 

2001 

Krause & 


when the bow shock leaves the grid sideways. Saxton et al. 
( 2002 |) have presented simulations with a closed injection 


velocity at kpc scale is not so easy measurable. Carilli & 


ratio. However, it should be emphasised, that this number 
is quite uncertain. 

Parma et al. (1999) measure the radio source life¬ 


Those studies have been reviewed 


by Normau (1993) and Ferrari (1998). One main result 
was that light jets, i.e. less dense than the surrounding 
medium, develop extensive cocoons, and are therefore be¬ 
lieved to be present in radio galaxies and quasars. The 
regime below that density, which will be called very light, 
was touched only occasionally. 

Recently, parameter studies of jets were carried out by 
Carvalho & O’Dea ( ^002a ,b|) and Saxton et al. ( 20021 ). 


These studies partly employ an open boundary condi¬ 
tion on the jet injection side. In that respect, they are 
complementary to the study presented here. The problem 
with the open boundary condition is that the backflow 
can transport large amounts of energy off the grid. This 
is a problem in the case of very light jets (especially for 
young bipolar sources), since the backflow is very fast com¬ 
pared to the head propagation. A similar problem arises 


time by synchrotron aging. Their straight forward ana¬ 
lysis yields a lifetime tg < 100 Myr for low power sources 
and tg < 10 Myr for high power ones. However, they point 
out that this discrepancy does not have to be real, since 
the estimate depends on the magnetic field, which is more 
uncertain in high power radio galaxies. Reducing the mag¬ 
netic field to one fourth of the applied equipartition value 
pushes the high power sources to 100 Myr, also. This 
would be supp orted by jet-to -counterjet length asymme¬ 
try constraints ( Scheuer 1991 ), and just be consistent with 
the age of the universe in the case of extended sources at 
the highest redshifts. Since there is the possibility that the 
older populations of synchrotron electrons are polluted by 
reaccelerated ones or mixed with younger populations in 
backflow regions, synchrotron ages indicate lower limits. 
Adopting 100 Myr as fiducial upper limit for radio galaxy 
lifetimes, Parma et alj (1999) arrive at a typical head ad¬ 
vance speed of Uhead = (0.5 — 5.0) kpc/Myr for the low 
power sources. It is not unreasonable that the same pa¬ 
rameters apply for high power ones. For those, ^cheuei 
(1995) gives an upper limit of Uhead = 30 kpc/Myr, from 
jet-to-counterjet length asymmetry measurements. 

The advance speed of the bow shock in jet direction 
can be computed from the one-dimensional momentum 
balance: 


side boundary and a big enough grid for the bow shock to 
remain inside for most of the simulation time. 

In this paper a systematic study of very light jets with 
77 G [10“^; 10“^] and internal Mach number Mj G [2.6; 26] 
is presented. The injection side boundary is closed, and 
the grid is big enough to follow the bow shock for the 
whole simulation time. Contrary to previous studies, the 
focus in this paper is on the jet physics specific for very 
low jet densities. First, observational constraints on the 
jet parameters are reviewed. It follows a description of 
the simulations. Afterwards, an analytic approximation 
for the bow shock propagation is given, and finally the 
results are discussed. The time unit used in the following 
is the Myr, an abbreviation for 10® years. 


2. Constraints from Observation 

In nearby radio galaxies, the movement of individual jet 
features can be observed in the inner 100 parsecs (e.g. 
Britzen et al, 2000] ). These measurements indicate appar¬ 
ent superluminal motion, which is a unique feature of rel¬ 
ativistic velocities, very close to the speed of light c. The 


Vb 


1 


1 + l/7in.^b 


V^V l-Fl/7jMj 


-f 1. 


( 1 ) 


Here, e is the ratio of beam to head cross section (simula¬ 
tions determine this parameter to e S [ 0 . 1 ; 1 ] for 77 > 0.01 
(e.g. Norman et al, 1983| )), Vb and Mb are velocity and 
Mach number (with respect to external medium) of the 
bow shock, e is a measure for the propagation efficiency. 
7 ni and 7 j are the adiabatic indices of external medium 
and jet beam, respectively. Equation (|^) can be approxi¬ 
mated for high Mach numbers and low jet density by 


77head = y/m Vj 


( 2 ) 


Barthel (1996) cite a value of uj « 0.4 c for the jet velocity 
in the kpc scale jet of Cygnus A, estimated from the hot 
spot spectra and consistent with luminosity, minimum en¬ 
ergy hotspot pressure, lack of internal Faraday dispersion 
in the lobes, and the jet-to-counterjet surface brightness 


77 is constrained by the fact that radio lobes develop and 
the above considerations (Eq. (||)): 3 x 10“® < 77 < 1. 
However, it will be shown below that Eq. (^ is no longer 
valid at these low jet densities. A lower limit of 10“^ for 
the density contrast is derived in Sect. ^ Unfortunately, 
the temperature in the jet beams is not known. Therefore, 
Mach numbers are unconstrained, observationally. 

Jet radii ( r,) are of the order of kpc (e.g. [Carilli fc 
Barthel 1996] ; Scheck et al. p002 ). From the total radio 
luminosities, one can estimate the total kinetic jet lumi¬ 
nosity. Since the source needs at least the energy to inflate 
the radio lobes, the total power has to be at least about 
twice the radio power (for Cygnus A, Carilli fc Barthe] 
1996), probably more. The most powerful sources at high 
redshift should therefore have Lkin > 10 "^® erg/sec. 
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3. Simulation Setup 


For an investigation of the parameter space of very light 
jets, axisymmetric hydrodynamic simulations were per¬ 


formed. The Newtonian 3D MHD code NIRVANA (Ziegler 


fe Yorke 1997) was used for the computations. NIRVANA 
is second order accurate and employs a monotonic upwind 
differencing scheme. It treats the following standard set of 
hydrodynamic equations: 


| + V.(pv) = 0 

( 3 ) 

— + V • (pvv) = -Vp 

( 4 ) 

de 

— -b V • (ev) = -p V • V, 

( 5 ) 


where p denotes the density, e internal energy density, v 
velocity, and p = (7 — l)e the pressure. Here 7 = 5/3 for 
a nonrelativistic gas is assumed. 

For a first scan of the parameter space 12 axisymmetric 
simulations were performed. The jets were injected in pres¬ 
sure equilibrium into a homogeneous external medium. 


Table 1. Simulation parameters 


V 

Vj [kpc/Myr] 

Mint 

-^ext 

10"^ 

95 

25.56 

255.6 

10"^ 

31 

8.083 

80.83 

10"^ 

9.5 

2.556 

25.56 

10"® 

95 

25,56 

808.3 

10"® 

31 

8.083 

255.6 

10"® 

9.5 

2.556 

80.83 

10~'‘ 

95 

25.56 

2556 


31 

8.083 

808.3 


9.5 

2.556 

255.6 

10"® 

95 

25.56 

8083 

10"® 

31 

8.083 

2556 

10"® 

9.5 

2.556 

808.3 


The hydrodynamic simulations are fully determined by 
the internal Mach number Mint, the density ratio p, and 
the pressure ratio, which is unity. Hence, they are scalable 
to the parameter range needed in a specific source. Table ^ 
gives Mach number (M), density contrast (p), and the ap¬ 
plied velocity. The density in the external medium was set 
to 2 cm“^, and the jet radius to rj = 1 kpc0. The grid of 
[600 X 600] points covered an area of Z x i? = 30 x 30 kpc^, 
so the jet radius corresponds to 20 grid points. Boundary 
conditions were set to axial symmetry on the axis, reflect¬ 
ing on the left-hand side, and open on the other sides. 
The reflecting boundary on the left-hand side is essential. 
If there was an open boundary, a substantial amount of 


The scaling was chosen in order to be easily comparable 
to previous work, which was a dapted to pa rameters assumed 
to be present at high redshift ( Krause| 2002). 


energy would leave the grid on that side. The simulation 
was stopped when one of the following events happened: 

1. The jet propagated to the right-hand side of the grid. 

2. The computation time exceeded a reasonable amount 
without promising new behaviour in the near future. 
Typically, this time was several weeks. 

3. The jet stopped at the nozzle, producing a shock at the 
inlet which ignored the shock jump conditions. This 
was caused by entrained material from the shocked 
external medium approaching the nozzle. 


4. Results 


4.1. Cocoons 

Contour plots of the final timestep for each simulation 
are shown in Fig. ^ As expected from higher ij simu¬ 
lations, the cocoons broaden with lower rj. Begelman & 
[Ciofli (1989) derive an analytic formula for the cocoon 
width rcoc- 


rcoc/rj « 


( 6 ) 


Figure shows that only the (M, 77 ) = (2.6,10“^) cocoon 
can be described properly by the above formula. For the 
other jets, the bow shocks are still too small to contain 
cocoons of the predicted size. Also, the geometry of the 
cocoon is quite different from the cylindrical one known 
from higher 77 jets. The density distributions show more or 
less spherical cocoons. Between 77 = 10“^ and 77 = 10“^, 
the cocoon undergoes a transition. Sometimes, the vor¬ 
tices it dissolves in are stringed in a line around the beam, 
but sometimes they join together forming a big vortex 
which extends approximately over the same size as the 
jet. Figure]^ (Af, 77 ) = (2.6,10“^) shows a new vortex just 
before being swallowed by the big one. A vortex with sub¬ 
structure is present in (M, 77 ) = (8,10“^). The interface 
between cocoon and shocked external medium suffers from 
Kelvin-Helmholtz (KH) instabilities. They start at small 
size in the vicinity of the jet head. As they grow, extend¬ 
ing their Angers into the cocoon, the backflow accelerates 
them towards the center. This behaviour is known from 


higher 77 simulations ( Krause fc Camenzind 2001 ) at high 
resolution. On the left-hand boundary, the Angers can get 
long enough to interact with the beam directly. In some 
cases, this caused a strong shock at the jet inflow, effi¬ 
ciently pushing the jet out of the computational domain 
(one of the reasons, why the simulation had to be stopped, 
see above). Quite often, this gas is drawn towards the jet 
head in between cocoon and jet beam, thus separating the 
cocoon from the beam. In the 77 = 10“® simulations, the 
initial conditions are not yet relaxed enough to determine 
the cocoon morphology reliably. 


4.2. Bow shocks 

The bow shocks in the presented simulations are quite dif¬ 
ferent from one another and from simulations with higher 
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Fig. 1. Overview over the 12 simulations performed in order to scan the parameter space. The Mach number (M) varies from 
2.6 to 26 from left to right. Density contrast ranges from 10“^ (top) to 10“® (bottom). The logarithmic density contours vary 
from dark (low) to white (high). The white spot in the center of the big vortex in the (M, 77 ) = ( 8 ,10“®) plot as well as several 
jet inlets are artificially white because of a problem in the plotting routine. In fact, they have the lowest values. The simulation 
time is indicated on top of the individual plots. 
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Fig. 2. The propagation efficiency measure t on average (left panel) and for the last 5% of the simulation time (right panel). 


77 . Morphologically, they can be classified in three groups. 
The extremes are located in the upper corners and the 
bottom row of Fig. |^. The jet in the upper right corner 
(M, 77 ) = (26,10“^) shows an elongated bow shock, com¬ 
parable to the simulation of Loken et al. ( |1992 ). Due to 
the high Mach number, the compression ratio is four over 
the whole surface of the bow shock. Reducing the Mach 
number to (M, 77 ) = ( 2 . 6 , 10 “^), the bow shock gets very 
regular (upper left corner). Evidently, this is caused by 
the lag of the contact discontinuity behind the bow shock. 
Even though there are some disturbances in the shocked 
external medium, they are not fast enough to catch the 
leading edge. The bow shock is still stronger in the di¬ 
rection of the jet propagation which produces a higher 
compression ratio near the head. 

If one would reduce the Mach number below one, no 
jet would form at all, and a spherical sound wave would 
propagate outwards. In the bottom row (77 = 10“®), the 
bow shock is spherical, reminding one of a spherical wave. 
Nethertheless, one can clearly tell from the density jumps 
that it is truly a bow shock. So, why is it spherical? A 
hint comes from the pressure versus number density (pn) 
histogram (Fig. From their start position at the ex¬ 
ternal pressure and the density of the jet and external 
medium ( 10 ® times jet density), respectively the fluid ele¬ 
ments move upwards to a thin line of high constant pres¬ 
sure. (For the (M, 77 ) = (2.6,10“®) and (M, 77 ) = ( 8 ,10“®) 
plots, the pressure increases so rapidly in the jet that the 
counts are too few to show up in the plots.) The pressure 
equilibrium of all the fluid parts affected by the jet can be 
understood, if one considers that the sound speed in the 
jet (roughly Ujet/ 2 ) is much higher than the bow shock 
propagation speed. Therefore, sound waves communicate 
and equalise pressure differences efficiently. 

The internal energy within the bubble of the bow shock 
is sufficient to drive a blastwave. The solution for the blast- 
wave is derived in Sect. ||. In this limit the jet propagates 
slowly, and mainly transforms its constantly provided ki¬ 
netic energy into heat. The resulting pressure drives the 
spherical blastwave. The velocity of the blastwave is pro¬ 


portional to . Hence, the blastwave is faster than the 
jet bow shock at early times, and slower later on. Equation 
( p^ determines the minimum bow shock radius at which 
deviation from the spherical blastwave can be expected. 
For 77 = 10“® it turns out to be 9.25 rj. The 77 = 10“® jets 
are already slightly above that limit. But the jets are in a 
quite complicated phase at the time shown in Fig. ^ For 
example in the case of (M, 77 ) = (26,10“®), one shock in 
the middle of the jet just became so strong that it devel¬ 
oped a backflow there. The old jet head at Z « 9 kpc is 
dissolving. This is a short phase, and because the blast- 
wave swept much of the matter away, the new jet will 
soon reach the position of the old one. Then, it will pre¬ 
sumably influence the bow shock. That this happens soon 
after the minimum radius of influence (ri) is reached, is 
evident from the 77 = 10 “'’’ density plots, where ri = Srj, 
and from the aspect ratios (see below). One can under¬ 
stand the bow shock shapes of the other simulations as 
combinations of these three types. 

The propagation efficiency measure e was determined 
from the simulations on average and for the final 5 % of 
simulation time. The result is shown in Fig. ||. The high¬ 
est values of e are reached for 77 = 10“®. In that case the 
bow shock is on average three times faster than the maxi¬ 
mum expected from the one-dimensional estimate. Again, 
this shows the high velocity of the blastwave. Also for the 
77 = 10 “^ and the 77 = 10 “® simulations the average ef¬ 
ficiency is still high, for the same reason. The 77 = 10“^ 
simulations have already reached an extension of roughly 
15 times their minimum interaction distance. Therefore e 
has fallen considerably below one. Also, it does not change 
significantly for the final 5 % of the simulation time. Slight 
variations in the propagation speed are typical for such 
jets. Extrapolating from higher 77 , the accuracy of e is a 


few p ercent for the given resolution (Krause & Camenzind 
^001 ). Summarising, very light jets blow up a big bubble, 
rapidly. After reaching ri, they stroll along with nonspec¬ 
tacular e. 

The aspect ratio of the bow shock was measured in 
regular intervals, and is shown in Fig. ^ as a function of 
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the position of the bow shock on the Z-axis. A general 
trend is evident for all simulations: Up to a propagation 
distance of two to three jet radii above r\ the aspect is 
approximately one. This is also consistent with the r] = 
10“® simulations. Then the aspect starts to rise. This can 
be explained by considering the bow shock expansion in 
the axial and radial directions. Up to ri, the bow shock has 
to expand according to the blastwave law, in all directions. 
After that phase, it changes to a typical jet propagation 
law, for the axial direction (compare Fig. |^). The jet head 
sometimes accelerates temporarily, which is well-known in 
jet physics {beam pumping, compare e.g. |Kossl fc Muller 
(|1988|)), but has on average a constant velocity, given by 

eVI). 

Concerning the radial direction, the bow shock is still 
pressure driven. It will be shown in Sect. I (Eq. (|^) that 
the Mach number of the bow shock at ri with respect to 
the external medium is given by the jet’s Mach number 
Mint- Hence, the bow shock is still strong at r = ri, at least 
for the M > 3 simulations. The deceleration is propor¬ 
tional to the aspect ratio via the pressure (assuming ellip¬ 
tical deformation), and should therefore be higher. On the 
other hand, secondary shocks in the shocked ambient gas 
can now reach the surface of the bow shock, accelerating it 
in limited regions. As an example, for the (M, ri)={8, 10 “^) 
simulation, the relative decrease of the radial bow shock 
velocity (1 — Vh/vg) between radial bow shock positions 
Va = 6.65 kpc and Vb = 17.0 kpc amounts to 47.8 %. 
One can easily derive from Eq. (0 that the blastwave 
law predicts a decrease of only 1 — (ra/ru)^/^ = 46.5 %. 
The agreement between the numbers is quite good in this 
case. In the other simulations there is more accumula¬ 
tion of jet material at the left boundary (compare Fig. |^, 
which sometimes strongly disturbs the bow shock near the 
left boundary. The difference between the measured and 
predicted velocity decrease can reach a factor of two there. 
However, the boundary condition at the left side compli¬ 



Time [Myr] 


Fig. 3. Bow shock propagation over time for the (M, rf) = 
(26,10“^) simulation. Shown is the position reached on the R 
and Z-axis, respectively. Logarithmic units were used for both 
axes. 


cates the interpretation. A preliminary conclusion is there¬ 
fore that the radial bow shock expansion keeps following 
the blastwave law for the simulated jets with M > 3 and 
r] < 10 “^ for the simulated times, if the disturbance by 
the material accumulated at the left boundary is not too 
strong. In a follow-up paper, where a 3D simulation with 
bipolar jets will be presented, this issue will be addressed 
in more detail. 

Consequently, the aspects increase with time. In Fig. |], 
the r] = 10 “^ and rj = 10 “"^ simulations have left the blast- 
wave phase not long ago, whereas the rj = 10 “^ jets are 
already more than ten times bigger than ri. The former 
jets are observed to push comparatively gently, basically 
elongating the bubble. The p = 10“^ jets show an interest¬ 
ing dependence on the Mach number: Only the low Mach 
number jet keeps the elongated spherical shape. 

The cocoons of the high Mach number jets are un¬ 
derexpanded according to Eq. (|^), and act violently on 
the bow shocks. In the phase where the bow shock is 
only slightly elongated, the bow shock velocity in axial 
direction has not yet reached the constant velocity. This 
can be seen from Fig. The gradual change between the 
two propagation laws is finished at t « 2 Myr, for the 
{M,r]) = (26,10“^) jet. At that time the bow shock has 
reached approximately 10 ri. The limiting aspect ratio (a) 
can then easily be computed from Eq. (H): 

a = = Mextv^i (7) 

Cs, ext 

where the index ext denotes quantities in the external 
medium. With the values from Table the limiting as¬ 
pect ratios turn out to be between 1.6 (M = 2 . 6 ) and 16 
{M = 26). None of the simulations has reached this limit 
yet. 


4.3. Beams 


A prominent feature of the jet beams is the strong shock 
roughly one jet radius behind the jet inlet (compare 
Fig. |). This indicates the first reflection point of the 
oblique shocks in the jet. The shock is caused by the high 
pressure and velocities of the backflow, which acts on the 
beam. The position of the first shock reflection point varies 
with time. This is consistent with simulation results at 
higher density contrast wit h closed left boundary condi¬ 
tions ( Kossl fc Miillei 1988 ). 

High Mach numbers (M > ^ can be sustained for the 
r] = 10“^ jets (compare Fig. With the exception of 
the {M,r]) = (26,10“^) jet, all the other beams struggle 
around M = 1 — 2. Quite often they become subsonic, even 
if they start with high Mach numbers. Strong shocks de¬ 
celerate those jets. Since also shocked ambient gas reaches 
the beam quite often, nothing prevents the beams from 
being disrupted by the KH instability. A fine example is 
the (M, 77 ) = (26,10“^) simulation (Fig. |^. The shocks 
in the beams sometimes get strong enough to shed vor¬ 
tices and produce backflows, thereby establishing a new 
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bow shock on Z axis [kpc] 



Fig. 4. Aspect ratio of the bow shock over propagation dis¬ 
tance. The upper panel shows the 77 = 10”^ simulations, the 
middle one the rj = 10 “® simulations, and the 77 = 10 ““* simu¬ 
lations are shown in the bottom panel. Logarithmic units were 
used for both axes. 


jet head. This behaviour is similar to axisymmetric sim¬ 


ulations of higher 77 jets (e.g. Kossl & Miiller 1988; Jones 


et al. 1399), but seems to be more violent at lower density 
contrast. 

The 77 = 10“^ jets suffer from the short evolution time. 
At the time shown, the jets still have the prominent shock 
at the jet inlet, which was used as initial condition. Since 
those are the computationally most expensive simulations. 


the problem can only be solved with a faster computer, 
evolving the jets longer. 

4.4. Pressure versus number density histograms 

From Fig. ||, it is evident that very light jets provide their 
own pressure everywhere. External medium as well as jet 
plasma get strongly shocked, and the high sound speed 
equalises the pressure within most of the jet bubble, espe¬ 
cially in the higher Mach number cases. The external pres¬ 
sure is therefore unimportant, and the pressure matching, 
used as initial condition, can be dropped in future simu¬ 
lations, without much difference. The pn-histograms show 
many counts at intermediate densities. They are caused 
by fluid elements resulting from mixing of jet plasma with 
shocked external medium. This is present in every sim¬ 
ulation, although there are typically less counts than in 
the region of shocked external medium densities (right) 
or jet densities (left). Especially prominent are several 
straight lines in the pn-histograms. They extend to very 
low densities, present in the center of vortices. The fol¬ 
lowing example examines the (M, 77 ) = ( 2 . 6 , 10 “^) jet. 
The density contours show two prominent vortices (dark) 
in the cocoon. Correspondingly, the pn-histogram shows 
two spikes extending down to a tenth of the jet density. 
Now, two areas, each centered on one vortex were cut out 
and examined separately. The individual pn-histograms 
are shown in Fig. Clearly, in the right figure, the up¬ 
per spike has disappeared. Nevertheless, it is prominent 
in the left figure. Hence, it is evident that vortices show 
up in the pn-histograms as straight lines. Therefore they 
follow the relation p (x nF. The left and right vortices of 
Fig. 0 have F = 0.35 and F = 1.26, respectively. The spikes 
in the other histograms have a maximum F of 1.5 — 1.7. 
If a vortex would originate from jet material with uni¬ 
form density and would be compressed in a shock with 
uniform compression ratio k, it would have a uniform en¬ 
tropy index s = p/rp. The pressure differences are caused 
by adiabatic changes due to the centrifugal force. The con¬ 
sequence is a p = srp subsystem, which is often observed 
in the pn histograms. The low F of the two vortices in the 
(M, 77 ) = (2.6,10“^) is exceptional. Systems of that kind 
are occasionally created in a strong shock in the vicinity of 
the jet axis connected with a temporal on-axis backflow. 
Such shocks are usually weaker in a 3D simulation. The 
life time of these systems is short. They mix rapidly with 
other material due to the limited grid resolution. 

The bow shock obeys a curved line in the pn- 
histogram. This is especially evident from the (M, 77 ) = 
( 8 ,10“®) simulation. The density contours show most 
of the space occupied by shocked external medium. 
In the pn-histogram, these fluid elements, starting at 
(log( 7 T,), log(p)) « (5,0), first move upwards, bending to 
the right. The maximum pressure and density in the 
shock is reached in a rightmost pinnacle ((log(7T,), log(p)) « 
(5.6, 2.2)). Then pressure and density decline fast to the 
isobaric regime. Here the density declines further. Even in 
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Fig. 5. Pressure over number density (pn) histograms. 100 bins were used on each axis. Pressure, number density, and counts 
are given in logarithmic units. Darker regions have more counts. Pressure and density have been normalised to the initial jet 
values. 
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Fig. 6. On-axis Mach number. 
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Fig. 7. pn-histograms for regions extracted around the two vortices of the (M, rj) = (2.6,10 jet. The left figure corresponds 
to the left vortex. 


this case, with the smallest cocoon, 87 % of the mass is 
concentrated in the region 8.5 kpc < r < 10.5 kpc, forming 
a thin shell. 

5. Analytic Approximation for the Blastwave 
Phase 

In the simulations of the last section the bow shock was 
found to resemble a spherical blastwave up to a certain 
time. The simulations show that the cocoons displace most 
of the external medium into a shell. Assuming a constant 
velocity v for a shell of mass M, one can find an analytic 
approximation for the expansion of this shell by postulat¬ 
ing force balance at the bow shock surface: 

^ {Mv) = S (Pint - Text) ■ (8) 

Here, S = 47rr^, and pint and Pext denote internal and ex¬ 
ternal pressure, respectively. The internal pressure is given 
by Pint = |(£'(i) — AIt'^/2)/(47rr^/3), where E(t) is the 
energy injected into the cocoon by the jet. 

Usually one considers self-similar solutions of the 
spherical force balance (e.g. Heinz et al. 199^ ; poker etlT 
2002; Krause 2002). Here, an alternative approach is pre¬ 


sented, computing the global integral of the equation. 

5.1. General solution 

In the case of negligible external pressure, which is the 
case for strong bow shocks, Eq. (0) can be integrated: 


j-r i-t 

/ M.{r)rdr = 2 / dti E{t 2 )dt 2 . 
Iq Jo Jo 


(9) 


density distribution, as long as both are integrable func¬ 
tions. 


5.2. Power law solutions 

A very useful type of solutions is found, if one assumes 
energy and density to be given by power laws: 

E = 

P = Po{r/roT 

The solution for the shell radius is in that case: 


_ / (k -I- Z){k + 5) r^C \ 
id -\- l)(d -|- 2) 27rpo/ 




( 10 ) 


For example, setting total energy and background density 
to a constant [k = d = Q) gives the well-known solu¬ 
tion for expanding supernova bubbles ( pedov 19591) . Jets 
can be assumed to deliver a constant amount of energy to 
their cocoons. Therefore d = 1 is appropriate. With that 
choice, the funct ional dependence of the self-similar ex¬ 
pansion law from Heinz et al. ( 1998 ) is recovered. However, 
for density exponents between -2 and 0, the constant of 
proportionality is by a factor of 4 to 1 lower in the present 
solution. This reflects the fact that it is, from the math¬ 
ematical point of view, a global solution. Therefore, for 
density exponents lower than -3 the bubble has to fight 
against an infinite amount of mass at r = 0. Equation (|^ 
consequently states that the blastwave will not expand at 
all. This is of course not found in a self-similar solution, 
since such a solution is always assumed to be far from the 
boundaries. A relativistic gas in the cocoon, as taken into 


Equation (||) is the general solution for any sort of energy 
input driving a blastwave into a medium with arbitrary 


account by Heinz et al. (1998), increases the constant of 
proportionality by roughly 10 %. Because the blastwave 
is decelerating for k > —2, whereas the bow shock from 
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a jet should propagate at approximately constant speed 
(compare Eq. (||)), the bow shock is first spherical, and 
gets shaped by the jet only if the jet thrust can push it 
faster. The expansion velocity of the blastwave v from a 
non-relativistic twin jet is given by: 


9k + 3 f ^ \ f ^ 
Ujy 4 {k + 5)^^° \ro ) 


( 11 ) 


where tjq = p-Jpa- Combining this with Eq. (^, and choos¬ 
ing To = J’j, results in: 


ri 


/16 (K-h5)4 
\8l {n + Sf 



O.5(?7oe")^ 


( 12 ) 


ri denotes the radius at which the jet first affects the bow 
shock. The latter approximation is valid for 0 > k > —2.9 
with an accuracy of 25 %. The coefficient goes to zero 
for K towards -3. The minimum ri is reached for e = 1. 
It may be surprising at the first glance that the jet head 
can even more easily reach accelerating bubbles. This is 
due to the fact that the jet head accelerates always faster 
than the blastwave, as long as the assumptions leading 
to Eq. (H) are appropriate. It follows that spherical bow 
shocks are expected for very light jets, in small systems, 
and in environments where the density does not decrease 
faster than with k « —2.9. The latter restriction does 
not apply, if the density exponent turns smoothly into -3, 
which is demonstrated below. 

The assumption of negligible external pressure holds 
as long as the bow shock is much faster than the sound 
speed in the external medium, which can also be shown 
explicitely from the equations in this section. The bow 
shock has a Mach number Mb with respect to the external 
medium at a radius r, in the isothermal case given by 



Fig. 8. Propagation of a blastwave powered by a jet into a 
King type atmosphere with (5 as indicated above. 


power law. A King type density distribution is often used 
for density distributions around galaxies: 


Pir) = Po 1 + 



- 3 / 3/2 


(15) 


The integrals in Eq. (|^) can be performed analytically for 
some cases of /3, resulting in: 


t = 


3 127rpo7’o 


F, 


(16) 


where Y is given by 

|a;(A^ -I- A/2) — iarcsinh(x)i3 
Y = ^ ^ (IxC — arctan(x)A^) 


iarcsinh(a;)C' — jxA 


> for 3/3 = 


-(17) 


with 


r 

ri 


/9 K + 3 -2/k\ 

V4(«:-b5)2'^° ) 




(13) 


Eor the conditions of the simulations presented here (k = 
0), Eq. states that the Mach number of the bow shock 
exceeds the jet Mach number when the bow shock reaches 
ri. If K < 0, the bow shock Mach number at r = ri is 
even higher for rjo < 10“^. It can reach typically several 
times the jet Mach number. In an isobaric atmosphere, 
the relation becomes: 



(14) 


where 6 is the overpressure ratio of the jet. Again, the bow 
shock is supersonic at ri for supersonic jets, unless they 
are heavily underpressured. 


5.3. King type solutions 


A = \/l -b a;2, 

B = 3/4-bx^ 

C = 3/2 -b x^, and 

X = r/ro. 


The solutions from Eq. ( |l7| ) are plotted in Eig. It is 
evident that the effect of a King atmosphere, which is a 
smooth transformation from a constant density profile to a 
decline with exponent —3/3, transforms the corresponding 
power law solutions smoothly into one another. Eor /3 = 
2/3, which corresponds to a /t = —2 power law solution for 
large x, the shell reaches constant velocity at infinity. For 
steeper density declines, the shell accelerates. In that case 
it is Rayleigh-Taylor unstable, and considerable mixing 
with the cocoon can be expected, as already pointed out 
by Heinz et al. ( 1998| ). In the same manner as for the 
power law case, one can compute ri for the King density 
distributions. For example, in the /3 = 1 case ri is given 
by the following formula: 


Extragalactic jets often propagate into environments 
where the density cannot be approximated by a single 




D{x,)-\ 


(18) 
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where 

xi = ri/ro, 

D(xi) = E(xi)~'^''^F(xi), 

3 

E(xi) = arcsinh(a:i)C'(xi) — —xiA(xi), and 

E(xi) = a:iarcsinh(a:i) — xf/A(xi). 

ri is plotted for three different core radii in Fig. ^ This 



Fig. 9. Minimum radius of jet influence (ri/rj) (horizontal 
axis) as a function of density contrast 77 = pj/po for different 
core radii ro in a /3 = 1 King profile. For ro = 00 the density 
is constant. 

density profile tends towards r~^ for large r. As predicted 
above, ri is finite also in that density regime. For rj > 10“"^ 
and ro > lOrj, the result is not much different from the 
constant density case. But m the more extreme cases ri 
can be reduced considerably, e.g. for 77 = iO ‘",ro = 2 rj, 
ri is reduced from 16 to 6 jet radii. 

Concerning the propagation law, the constant density 
formula approximates the King solution to 10 % up to 
approximately one core radius for the shallowest density 
profile and half a core radius for the steepest one (Fig. 

Summarising, the constant density formulas can be 
used to approximate the true solutions if ri is not much 
bigger than the core radius, if the density contrast is above 
10 “^, and the density profile steepens not much more than 
into a power law. 

6. Discussion 

A very interesting feature of the above results is the blast- 
wave phase. Since the blastwave is faster at early times, 
it paves the way for the jet, effectively increasing the jet 
propagation velocity. Later, the jet head is faster than the 
blastwave equation predicts. Therefore, an upper limit for 
jet ages is given by Eq. This can be used to estab¬ 
lish a lower limit for 77 . For very light jets, one can regard 
ri as a typical extention. The time for a jet to reach the 
distance ri in a constant density atmosphere is: 


t{ri) = 


1 


2^ 




(19) 


With the observational constraints cited in Sect, g, one 
obtains a new lower limit for the density contrast in ex- 
tragalactic jets: 77 > 10 “”^. 

It was shown that the aspect ratio of the bow shock 
depends on the source size. It is close to unity in the blast- 
wave dominated phase, and starts to increase soon after 
a critical size (ri) is reached. This increase will come to 
an end when the sideways expansion velocity drops to the 
sound speed, ^axton et al. ( 2002 ) have computed simi¬ 
lar models with density contrasts down to rj = 10“^. The 
simulation results are generally comparable in the regions 
where the simulation parameters overlap. The aspect ra¬ 
tio at the end of their 77 = 10 “^ simulation is also close 
to 1.2. The aspect ratio can be used to constrain density 
ratios. For example, the aspect ratio o f 1.2 for the bow 
shock in Cygnus A ( Smith et al. 2002|) together wit h its 
extension of roughly 120 Cj ( Carilli fc Barthel 1996|) are 
clearly inconsistent with 77 > 10“^. Probably, 77 should be 
below 10“^. With cluster temperatures of approximately 
10® K ( ^mith et al. 2002 ), resulting in probable exter¬ 
nal Mach numbers of a few hundred, one derives a limit¬ 
ing aspect ratio of about 2 , which the source has not yet 
reached. However, the jet clearly acts on the bow shock, 
and therefore one can check Eq. dl) , which is fulfilled if 
77 > 3 X 10“®. Another result was that low cocoon width 
to bow shock width ratios can only be achieved in low in¬ 
ternal Mach number jets. This result is depends crucially 
on the closed left boundary, which is essential in order to 
keep the internal energy inside of the computational do¬ 
main. Other studies have shown that the cocoon can be 
small even in high Mach number jets for an open bound¬ 
ary on the injection side ( [Carvalho fc O’Dea 2002a; Saxtou 
et al. 2002). Applying Eq. (|^ results in a probable range 
of 77 G [10“®; 10“®] for Al G [2; 6 ] for the jet in Cygnus A. 

With the general solution of the spherical force bal¬ 
ance, one can also predict how the results achieved here 
change in a nonhomogeneous density distribution. It was 
shown that - for a King profile - the blastwave phase is 
a little shorter for very small core radii. But in most real 
cases, the blastwave phase will remain up to similar source 
sizes as in the constant density case. This will be similar 
in all profiles with a constant density in the central region. 
Preliminary results from simulations with a King density 
distribution support this conclusion. 

The presented results are subject to two potentially 
important limitations: The axisymmetry and the left hand 
boundary. 3D simulations of light jets have shown that big 
vortices are unlikely to survive ( Tregillis et al. 2001 ). One 
should therefore expect turbulence on smaller scales in 
the 3D case. The jet beams should be less disturbed, and 
might be able to sustain higher Mach numbers. So, they 
might infiuence the bow shock earlier than the 2D jets. 
However, since the 2D jets show this interaction already 
shortly after reaching the theoretical ri, the difference 
should not be severe in that respect. The left boundary 
is an obstacle for material that would otherwise leave the 
grid there, and interact with the counterjet. This should 
change the appearance of the cocoon in regions around 
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Z = 0. Both these issues will be addressed in a future 
paper where a 3D simulation with bipolar (back-to-back) 
jets will be presented. 

By now, more and more X-ray bubbles are detected 


around radio sources, mainly in CD galaxies (e.g. Soker 
et al. 2002, and references therein). Such bubbles can be 


interpreted as gas contained within the bow shock (or the 
bow wave, if the shock has relaxed) from the jet. The the¬ 
ory presented in this work should be applicable to some of 
those sources. The typical bubble size for which one should 
expect supersonic bow shocks (homogeneous atmosphere) 
is given by: 


r = 24 


0.1cm ^ £ /lOOOkm/s^^ 

no 1046erg/s V Cs 


1/2 

kpc, (20) 


for typical values of the sound speed Cg and number den¬ 
sity no in the external medium. We should therefore ex¬ 
pect no strong bow shocks around jet sources with tens of 
kpc or more in diameter (compare e.g. in Hydra A, Nulserl 
et al. 2002). This conclusion is not changed if one consid¬ 


ers elongated bubbles, since the increased volume due to 
the elongation can only decrease the pressure inside the 
bubble. Hence, Eq. |2^ gives an upper limit. This means 
that the main energy transfer from radio sources to the 
gas of a galaxy cluster happens at early times in the jet 
evolution. The power input into the cluster gas is 30 % of 
the total jet power for the constant density case. This fol¬ 
lows from integration of the energy gain at the bow shock. 
A strong jet source therefore has delivered 


Egas = 10®°erg 


f. 


C 


3/2 


\ 10^® erg/s^ 

/ r] / »o \ 

Vl0“"^/ VO.lcm”®/ 


c/3 

- 1/2 


- 5/2 


( 21 ) 


to its surroundings when its bow shock reaches ri. This is 
already a considerable fraction of the total energy of the 
gas in a typical galaxy cluster. Cooling flows can therefore 
easily be stopped by a powerful jet. One such episode ev¬ 
ery 1000 Myr could balance the energy loss and keep the 
cluster temperature in the observed keV range. 

The results suggest that spherical bubbles should be 
found around most of the extragalactic jet sources at early 
evolutionary stages. The classes of compact symmetric ob¬ 
jects (CSO) and GHz peaked sources (GPS) are associated 
with such sources. Unfortunately, the bubbles are quite 
small. In order to resolve a region of a few kpc in diameter 
with Chandra, the source has to have a redshift consider¬ 
ably below one, where few objects where found ( pnellerj 
& Schhizzi 2002| ). The hot (10® K) bubble would radiate 
in the X-rays at about 10^^ erg/s, hard to discern from a 
quasar background. 

What about the radio morphology? As was shown 
above, young jets in the bubble phase are often disturbed 
or disrupted. This is also true for CSO/GPS sources 
( ^aikia et al. 2002). They are e.g. much more asymmetric 
than large scale jets. An old puzzle with CSO/GPS sources 


was if they were stuck within their galaxy because of an es¬ 
pecially dense interstellar medium (ISM) or because they 
are young. With the results obtained above, one can ex¬ 
clude the first possibility. Let us consider for example a jet 
with a power of 10^^ erg/s propagating at near light speed 
through an incredibly dense ISM, 10® times denser than 
the jet. For the first kpc, a typical core radius, it would 
take the jet head about 30 Myr, if propelled by the direct 
thrust alone. Such a jet was probably in the bubble phase 
for most, or maybe all, of its lifetime. The upper limit on 
its age from Eq. (0 turns out to be 0.3 Myr, even for 
an ISM density of 100 cm“®. Therefore, the probability of 
finding a source at a small angular size is not increased 
by considering dense environments. Recent age determi¬ 
nations by spectral aging methods and VLBI kinematics 
of hotspots str ongly argue for young source ages of CSOs 
( Conway 2002 ). 

The KH instability at the contact discontinuity is ubiq¬ 
uitous. This is especially interesting since the same be¬ 
haviour was found for ry = 10“^ at highest resolution 
(Krause & Camenzind 2001). The KH instability entrains 
shocked external medium into the jet cocoon, growing and 
accelerating it towards the center. This mechanism seems 
to be a promising candidate for bright X-ray features or 
the produc tion of emission line regions, where the densities 
are higher ( Bicknell et al. 200C): the gas is accelerated and 
pulled down inside the cocoon, where it can radiate by re¬ 
combination or reradiation of quasar light. Given the high 
energy content of the radio lobes, one could also imag¬ 
ine that they somehow transfer energy to the entrained 
gas. Since the cocoon plasma is quite exotic, the classical 
heat conduction formulas are probably not appropriate. 
Therefore, exact computations are beyond the scope of 
this work. 

Goncerning the beams, one can fairly say that those 
non-relativistic jets without significant magnetic field can¬ 
not travel to large distances without being disrupted. One 
needs a magnetic field in order to keep some of the co¬ 
coon around the beam, and to damp the KH instabilities 
in the beam. Alternatively, they would always look very 
disturbed, morphologically, until the bow shock had swept 
enough material away, so they could effectively travel into 
a less dense medium. Since bow shock detections point 
clearly in the direction of strong density contrast (i.e. 
77 ^ 1 ), a solution to this problem can be expected from 
further research. 
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